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We consider the pairing Hamiltonian and systematically construct its density functional in the 
strong-coupling limit and in the limit of large particle number. In the former limit, the functional 
is an expansion into central moments of occupation numbers. In the latter limit, the functional is 
known from BCS theory. Both functionals are nonlocal in structure, and the nonlocalities are in 
the form of simple products of local functionals. We also derive the relation between the occupation 
numbers and the Kohn-Sham density. 
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PACS numbers: 21.60.-n,21.60.Fw,71.15.Mb,74.20.-z 
I. INTRODUCTION 

Density-functional theory (DFT) is a very popular the- 
oretical method in many branches of physics, as it allows 
one to obtain the ground-state energy of an interacting 
many-fermion system from the solution of a mean-field 
equation 0, |2( • The main challenge in DFT is the con- 
nection between the Hamiltonian and the form of the 
density functional. Usually, one approximates the un- 
known (and supposedly nonlocal) density functional in 
terms of local densities, currents, and gradients of the 
density. This approach has led to impressive results in 
quantum chemistry and nuclear physics 0, HI IE • 

Pairing plays an important role in atomic nuclei (see, 
e.g., Ref. |8| for a recent review). Within the mean-field 
approach to nuclear structure, pairing is often included 
in the form of the BCS model added to the Skyrme func- 
tional 0, or th roug h Hartree-Fock-Bogoliubov (HFB) 
calculations [Tol lllj] The formal extension of Kohn- 
Sham DFT to superconductors also leads to HFB equa- 
tions [l^ ■ I n these latter two approaches, the energy is a 
functional of the density and the pairing field (anomalous 
density). It is thought that the latter is a highly nonlocal 
functional of the density and thereby captures aspects of 
the density functional that would be difficult to model 
otherwise. Within the HFB approach, the particle num- 
ber is not a conserved quantity anymore, and particle 
number projection becomes a concern and additional ex- 
pense 13]. In this work, we revisit the pairing problem 
and construct its density functional without resorting to 
the HFB approach. We thereby avoid some of the prob- 
lems associated with pairing fields and gain insight into 
the nonlocal structure of the density functional. 

The direct and systematic construction of a density 
functional from a given Hamiltonian is possible only for 
solvable or sufficiently simple systems. We mention, for 
instance, the systematic construction of the density func- 
tional for a dilute Fermi gas with repulsive interactions 
[Til Il5l ] , or the description of the Fermi gas close to the 
unitary limit in terms of local densities |l(j and gradient 
corrections 0,0]. Albeit being simple, these systems 
are nontrivial, and they provide us with useful insights 



for the empirical construction of the nuclear density func- 
tional. Examples are, e.g., the occurrence of terms pro- 
portional to p 7 / 3 0] or proportional to p 5 / 3 ^3] which 
are absent in presently employed functionals. The pairing 
Hamiltonian is another solvable problem, and the deriva- 
tion of its density functional is the purpose of this paper. 

This paper is organized as follows. In Sect. EI we de- 
scribe the pairing Hamiltonian and the method to derive 
its density functional. In Sect. IIIII we construct the den- 
sity functional in the limit of large numbers of pairs. In 
Sect. llVl we construct the density functional in the limit 
of large coupling strength. Both density functionals will 
be given in terms of occupation numbers. In Sect. we 
work out the relation between occupation numbers and 
the spatial density, and discuss extensions of presently 
employed density functionals to include pairing. We close 
with a summary of our results. 



II. HAMILTONIAN AND METHOD 

The pairing Hamiltonian is defined as 

n n 
H = Y,Y, £ J & ls & js a i; a iT 6 J'T a il • (!) 

J=l s=J.,t i,i=l 

We thus have fl single-particle orbitals which can be dou- 
bly occupied with a constant pairing interaction. This 
pairing model was studied and solved by Richardson in 
the 1960s as an alternative to the BCS approximation 
jl9( . Recently, the Richardson model has found renewed 
interest in the framework of superconductivity in small 
mesoscopic systems, and we refer the reader to the review 
by Dukelsky et al. [2(|. The Hamiltonian JIJ is exactly 
solvable and integrable (i.e., there are as many conserved 
quantities as degrees of freedom |2l|). The ground-state 
energy 



N 



(2) 
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of N pairs of fermions is given in terms of the solutions 
Ei of the Richardson equations 



N 

E 

3=1 



En — Ea 



fj 



-E 

fc=i 



i 



i=l,...,N . (3) 



Note that the ground-state energy J2J depends implicitly 
on the single-particle energies Note also that the 
ground state of the pairing model is a superposition of 
fully paired states, since singly occupied orbitals are not 
subject to the interaction. 

We want to construct the density functional for the 
pairing Hamiltonian. This requires us to compute the 
Legendre transform 



F (l n k}) = E- ^n k e k 



(4) 



k=l 



of the ground-state energy J5J) with respect to the occu- 
pation numbers 



dE 

de k 



(5) 



Due to the Hellman-Feynman theorem, the occupation 
numbers JSJ are indeed the expectation value n k = 
(V>l(SjL fife; + SjL a k ^)\ip) of the ground-state Usu- 
ally, DFT practitioners work with the density instead 
of occupation numbers. However, the Hamiltonian is 
given in terms of the Fock-space operators a\. and a k , and 
its density functional is therefore naturally expressed in 
terms of occupation. The relation between the density 
and the occupation numbers will be given in Sect. Ivl 

In order to actually perform the Legendre trans- 
form (0} , we need a closed expression of the ground-state 
energy in terms of the single-particle energies. Unfortu- 
nately, no such expression is available for the solutions 
of the Richardson equations , as these equations have 
to be solved numerically [22, |23|. However, closed ex- 
pressions for the ground-state energy do exist for the two 
limiting cases of a large number of pairs N S> 1 or a large 
coupling strength g/e k ^> 1. In Sect. 11111 and Sect. II VI 
we will construct the density functional of the pairing 
Hamiltonian in these two limits, respectively. 

Note that the density functional is only determined up 
to an overall constant. Shifting the single-particle ener- 
gies e k by a constant does not change the corresponding 
occupation numbers n k . Indeed, summing Eq. JSJ over 
k yields n k — 2N, which is independent of the 

single-particle energies e k . Thus, it is not possible to in- 
vert the equations 10 and to express the single-particle 
energies in terms of the occupation numbers. To avoid 
this problem, one can introduce an orthogonal transfor- 
mation from the single-particle energies e k to new vari- 
ables ei,...,en such that e\ = Y^ k =i £ k- The 
new variables ea, . . . , en are invariant under a constant 
shift of the single-particle energies e k (since they are or- 
thogonal to ei). The corresponding new occupation num- 
bers v k — dE / de k are related to the occupation numbers 



n k by the same orthogonal transformation. One has in 
particular V\ = 2N/Vfi and F({v k }) = E — J^k^k^k- 
The inversion can be performed for the new variables 
e2, ■ ■ ■ , en. Note that variations of the resulting density 
functional F({i> k }) with respect to U2,---,vq conserve 
the number of particles. 

For the cases considered in Sect. IIIII and in Sect. IIVI 
respectively, we can follow a simpler path. The ground- 
state energy depends only on certain combinations of 
single-particle energies. These can directly be expressed 
in terms of combinations of occupation numbers without 
solving the individual e k in terms of the n k . 



III. DENSITY FUNCTIONAL IN THE LIMIT 
OF LARGE PARTICLE NUMBER 

In the limit of a large number of pairs N ^> 1, the 
Richardson equations can be solved in a systematic 
expansion 0,E3- This leads to the BCS energy 



Ebcs = Y, £ k-^-2N) + ^( £ fe ~ + A ' ■ 

k=l " k=l 

(6) 

The gap A and the chemical potential fj, determined by 

n 



E 



i 



9 fr[ ^(e k - ^ + ' 



fi - 2N = 



£k - /-* 



(7) 



(8) 



We are interested in the construction of the density 
functional for the pairing Hamiltonian with ground-state 
energy © and perform the Legendre transform |Q in 
three steps. First, we compute the occupation numbers 



dE 



BCS 



3ea 



= 1 - 



(9) 



Here, we made use of Eqs. © and 0. The form Eq. © 
is intuitively clear: For vanishing gap, the occupation 
numbers are two and zero for orbitals below and above 
the chemical potential, respectively. Second, we invert 
this equation and obtain the single-particle energies as a 
function of the occupation numbers 



(1 



n,j)A 



y/ nj {2-n 3 ) 



(10) 



Here, the gap and the chemical potential are also func- 
tions of the occupation numbers. Indeed, summing 
Eq. 10 over all orbitals yields Eq. JSJ), while we obtain 



A = f 



(ii) 



3 



from the gap equation J7J. Third, we compute the Leg- 
endre transform and obtain the density functional 



u 

Fbcs = Ebcs — ^ < 

3=1 

A 2 




(12) 



The equations presented in this section are well known 
from BCS theory [9j, |26j. The form of the functional 
clearly exhibits two properties of the pairing interaction, 
(i) Pairing creates an instability of the Fermi surface, 
and (ii) is nonanalytical for weak couplings. Property 
(i) is included in the functional (|12[1 as there is a gain 
of energy associated with deviations of the occupation 
numbers from two and zero. Property (ii) is also fulfilled 
since the functional is nonanalytical for small values of 
2 — rik and rik, respectively. Note that the functional is 
nonlocal, as it is the square of a local functional. Note 
also that the functional Ijl2(l is among the simplest ones 
that includes nonlocality, overall quadratic scaling in oc- 
cupation numbers (motivated by a two-body interaction) , 
and properties (i) and (ii) . From this point of view, the 
form of the functional could almost have been guessed. 



IV. DENSITY FUNCTIONAL IN THE 
STRONG-COUPLING LIMIT 

Pairing in nuclei is often associated with the weak- 
coupling regime. An exception are, e.g., nuclei with 
magic proton number where the neutrons partly fill a 
single j-shell. For such nuclei, pairing is dominant due 
to the degeneracy of the single-particle orbitals. We con- 
sider the corresponding strong-coupling regime in this 
section and gain additional insight in the nonlocal struc- 
ture of density functional. 

Within the strong-coupling regime, the ground-state 
ener gy can be computed within the quasi-spin formal- 
ism [23 - Here, we follow approach by Yuzbashyan et 
al. |25j. In the strong-coupling limit, the pairing term 
is the dominant part of the Hamiltonian Q, and one 
has g ^ \sk — e|. Here, and in what follows, we denote 
averages over single-particle orbitals by the overbar, e.g. 



(13) 



k=l 



In the strong-coupling limit, the energies Ei fulfill Ei 3> 
and one can expand the Richardson equations @ 
in terms of the small parameters £k/Ei. This yields 
sums of the form . E~ p which in turn are expanded 
into a series of inverse powers of the coupling g. The 
corresponding expansion coefficients can be determined 



through recursion relations. We formally extend the re- 
sults by Yuzbashyan et al. |25| and write the the ground- 
state energy as 



E({e J })=E^+2Ne + 9 J2Y,^jI[( £ -^ ■ (W) 



CA TT / -nA 

j>2 A " <=1 

Here, N denotes the number of fermion pairs, and 

= -N(fl - N + l)g (15) 

is the leading contribution of order g to the ground-state 
energy. The sum over A in Eq. I|14[l runs over all par- 
titions A = [Ai, A2, . . . , Afe] of j into integers Ai > A2 > 
• • • > Afe > 2 such that 



5> 



(16) 



Let us analyze the ground-state energy 1)14(1. The first 
term is of order g and is quadratic in particle number. In 
the strong coupling limit, each single-particle orbital has 
approximately equal occupation 2N/CI, and the second 
term on the right hand side of Eq. i|14|) sums the corre- 
sponding single-particle energies. The subsequent terms 
contain ratios of single-particle energies and the pairing 
strength, and these central moments are invariant under 
shifts of the single-particle energies. 

Yuzbashyan et al. estimated that the expansion l|14|) 
converges for ratios (sj — e)/g ~ 1. They gave explicit 
expressions for the first two coefficients c\. Snyman and 
Geyer j2^ corrected a typo in the second coefficient, and 
Barbaro et al. |30| gave an expression for the two fourth- 
order terms. We confirmed these results by following the 
steps described in Ref . j2^| . The first few coefficients are 
thus 



C[2] = 

C[3] = 

C[4] = 

C[ 2 ,2] = 



(17) 
(18) 



4A(fi - N) 
1) ' 

8N(fl~ N)(n-2N) 

n 2 (n-i)(n-2) ' 

16AT(Q- N) 
~Q 3 (n- l) 2 (ft-2)(f2-3) 
x (n 2 (Q- 1) - N(Q- iV)(5fi-6)) , (19) 
16N(Q — N) 



n 3 (n-i) 2 (n-2)(n-3) 

-N(tl-N) (3(30-4) - 



n 2 {2fl -3) 
CI 



n- 1 



, (20) 



and higher-order coefficients can be worked out by fol- 
lowing the recursion relations given in Ref. |2fij . 

We are interested in the construction of the density 
functional Q for the pairing Hamiltonian with ground- 
state energy E. This requires us to perform a Legendre 
transform of the ground-state energy (|14|) with respect to 
the single-particle energies (i.e. the external potential). 



4 



This transformation is done in three steps. First, we 
compute the occupation numbers 



dE{{e k }) 



(21) 



This is most easily done by utilizing the identity 

^-(^^(fo-D'-'-^iP). (22) 

Second, we express the single-particle energies as func- 
tions of the occupation numbers. This inversion is most 
conveniently done order by order in terms of the strong 
coupling expansion. Third, we have to compute the sum 
on the right hand side of Eq. Q. To this purpose we 
employ the identity 



n 

E n i s i 

3 = 1 



n 

Ec 



i)(e j -e) + 2Ns , (23) 



where 



2N 



(24) 



is the average occupation number. 

Let us follow these steps order by order in the strong 
coupling expansion. In leading order (LO), the ground- 
energy is 



E w ({e 1 })=E^ +2Ne : 
and the occupation numbers (J3J) are 



(25) 



(26) 



This expression is independent of the single-particle en- 
ergies, and the Legendre transformation cannot be per- 
formed. The formal reason is, of course, that the en- 
ergy H25|) does only depend on the sum of the single- 
particle energies (cf. the discussion at the end of Sect. lTT|l . 
In the strong-coupling limit, all single-particle orbitals 
have equal occupation. Thus, the individual occupation 
numbers do not depend on the single-particle energies. 

In next-to- leading order (NLO), the ground-state en- 
ergy is 



C \2] 



E^o({e 3 }) = E w ({e 3 }) + -^{e-eY . (27) 



The occupation numbers fulfill 

- ^[2] 



(e 3 - e) 



(28) 



Note that cpi < 0, i.e. the occupation of orbitals be- 
low (above) the average single-particle energy is above 
(below) the average occupation number. For the inver- 
sion, we take the square of Eq. (|28(1 and average over the 
single-particle orbitals. This yields 



7 =72 gW~ =72 

(e-e) =- (n-n) 

4croi 



(29) 



Finally, we compute the sum 



E 



2c 



n,-e,- 



[2] 



.9 



(e - e) + 2Ar e 



(30) 



Putting all together, we arrive at the density functional 
in NLO 

o 

-FNLO ({"■&}) = ^NLO — ^ Hfcgfc 

fe=l 



i) g^ 2 



(n-n) 2 . (31) 



4c 



[2] 



In the presence of an external potential, given in terms 
of single-particle energies e 3 - , the density functional thus 
becomes 



-Enlo ({rifc}) = F NLO ({n fe }) 



n 

E< 

4=1 



(32) 



The ground-state energy for the system with N pairs 
is found from the requirement that number-conserving 
variations of the functional vanish. The corresponding 
occupation numbers fulfill 



dE Nho ({n k }) 
drii 







(33) 



and it is understood that the derivative does not act on 
the iV-dependent constants and cr* of the func- 

tional JSU- It is straightforward to check that the so- 
lution of this equation yields the energy l|27|) . 

For the calculation of the density functional in N 2 LO 
we start from the expression 



E 



:,,■:{ ;}) = £NLo(fe}) + ^(e - ef (34) 



for the energy, and repeat the three steps outlined above. 
For the inversion, we raise the expression corresponding 
to Eq. H28f) to the second and third power, and average 
over single-particle orbitals. We neglect any terms of or- 
der O ([(sj — s)/g] 4 ) or higher that are generated in this 
procedure. This yields two linear equations that express 
the central moments of occupation numbers (n — n) for 
fi = 2,3 in terms of the central moments of the energy 
(e — e) M . The inversion is straightforward, and the final 
result is 



F N 2 LO ({n k }) = F NLO ({n fc }) + 



gc[ 3 ] 



«4] 



(n-nf . (35) 



Let us pause for a moment and consider Eqs. I|31(l and 
iJSSJ, respectively. The strong-coupling expansion of the 
ground-state energy (f 141) . given in central moments of 
the single-particle energies, translates into an expansion 
of central moments of occupation numbers. Clearly, devi- 
ations of the occupation numbers from their mean value 



5 



n increase with increasing value of the ratios (ej — e)/g. 
In the limit of vanishing pairing strength g — > 0, the 
occupation numbers become discontinuous at the Fermi 
energy. Note that the density functional is local in the 
occupation numbers up to N 2 LO, as only simple sums 
over occupation numbers appear. Nonlocalities start to 
appear at N 3 LO, and we will address this order in what 
follows. 

At N 3 LO, the ground-state energy is 



up to the power of j, i.e 



E^Lo{{£j}) = 



E N2LO ({e 3 }) + ^(e-ef 



'[2,2] 



9° 



{e-ef 



(36) 



As we will see, the last term in this equation gives rise to 
the nonlocality in the density functional. Such products 
of central moments of single-particle energies are ubiqui- 
tous at higher orders. For the occupation numbers, we 
find 



2c 



n,- — n 



[2] 



+ 



gQ 
3c [3] 



(ej - e) 



g- 

4C[ 2 ,2] 



( Sj -e)(e-eY 



(37) 



We raise this equation to the power two, three, and 
four, average over single-particle orbitals, and neglect 
any terms on the right hand sides that are of order 
O — £)/g] 5 ) or higher. The resulting three equa- 
tions express the central moments of the occupation num- 
bers as nonlinear functions of central moments of the 
single-particle energies, and the nonlinearities are due to 
the terms proportional to croai. We invert these equa- 
tions and neglect terms that are of order O Urij — n) 5 ) 
or higher. The final result for the density functional is 



F({ nj }) = + gJ2J2 d ^ j U( n -^ Xi ■ ( 39 ) 

j>2 A i=l 

Here, A again denotes the sum over partitions (cf. 
Eq. (|16|) ). and coefficients d\ are functions of the co- 
efficients c\ in Eq. I|14l) . The functional is nonlocal, but 
the multiple sums over occupation numbers simply fac- 
tor into products of simple sums. Thus, the nonlocal 
structure is rather simple. It is also inexpensive from a 
computational point of view. 

It interesting to compare the structure of the func- 
tional l|39|) in the strong-coupling limit with the func- 
tional from BCS theory. Within the BCS approxi- 
mation, the nonlocality is limited to the square of a local 
functional. This suggests that products of multiple local 
functionals, as they appear in (|39|l . must be corrections 
of order 1/iV or smaller. This can be verified by expand- 
ing the BCS functional l|12[l around the strong coupling 
limit. We thus set nj = n + (rij — n) in Eq. l(T2")l and 
keep only terms up to second order in n.j — n. The re- 
sult agrees with the NLO result 1|31[) up to a term —gN. 
This term is recovered when corrections of order 1 /N to 
the energy © are included (see, e.g. the first term of 
Eq. (2.16) in Ref. [28|). 



V. FROM OCCUPATION NUMBERS TO THE 
DENSITY 

How can one extend presently employed density func- 
tionals 0, |32 to include pairing effects? A number of 
approaches can be found in the literature 

USUI. Here 

we follow the direct path that links the occupation num- 
bers to the density. Recall that DFT is built on the 
ground-state density 



p{x) 



«=bT 



(40) 



F N 3 LO ({n k }) = 



F^2 ho ({n k }) 

4C[ 2 ]Cp,2] 



yC [3] 



64cS 



4c 



[2]C[ 4 ] 



Hi 



64cf 2] 



gVl (n — n) 



(38) 



The second term on the right hand side of this equation 
is nonlocal in the occupation numbers, since it involves a 
double sum. This form of nonlocality is, however, rather 
simple, as the double sum is a product of two individual 
sums. 

We can now generalize the expansion of the density 
functional to higher orders. The density functional at 
order j consists of products of central moments of the 
occupation numbers, and includes occupation numbers 



Here, ^(x) creates a fermion with spin projection s at 
the position x, and denotes the ground state. How- 
ever, the energy functionals l|12|l and l|39|l are based on 
the occupation numbers n k that specify the occupation 
of the single-particle orbitals of the pairing Hamiltonian 
We denote the corresponding single-particle wave 
functions as uu{x) and have 



d 3 x u k (x)^l{x) 



(41) 



Within Kohn-Sham DFT, the ground-state den- 
sity (|40[l is given in terms of the occupation numbers 
n a of the Kohn-Sham orbitals <f) a (x), i.e. 



p(x) 



]n a (j)* a (x)4> a (x) 



(42) 
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Here, n a = 2 (n a — 0) for occupied (unoccupied) Kohn- 
Sham orbitals in fully paired systems, and we use the con- 
vention that Greek indices denote Kohn-Sham orbitals 
while Roman indices label the single-particle orbitals of 
the pairing Hamiltonian ifl}. The occupation numbers n k 
are related to the Kohn-Sham orbitals and Kohn-Sham 
occupation numbers via 

nk = }^{k\a)n a (a\k) , (43) 

a 

with 

(a\k) = J d 3 x<f>* a (x)u k (x) . (44) 

Thus, the occupation numbers l|43() are functionals of the 
Kohn-Sham orbitals and therefore nonlocal functionals of 
the density. The functional derivative of interest is 

SF({n k }) = y, dF Sn k 
6<fc(x) ^ dn k 6<j>* a (x) 

= n a }- — (k\a)u k (x) . (45) 
z r' on k 

k 

Thus, one might add a functional of the form 13911 or 112(1 
to commonly employed density functionals [3lL l32j . The 
single-particle wave functions u k {x) are determined by 
the pairing Hamiltonian; they can be, e.g., plane waves or 
shell model orbitals. The nonlocality of the density func- 
tional has two sources. The first is the nonlocal structure 
of the functionals (|39|l and ((12(1 in terms of the occupa- 
tion numbers, and the second is due to the relation i|43|) 
between the occupation numbers and the Kohn-Sham or- 
bitals. 

The use of occupation numbers (|43|l in a density func- 
tional is only new at first sight. The Kohn-Sham kinetic 
energy 

h 2 r 

T = ^$>a / <PxV4>* a {x)-V4> a {x) 

a 

= Y,^— n k 46 
^— ' 2m 

k 

can also be written in terms of occupation numbers n k 
of plane wave states. Equation 1(46(1 is obtained through 
integration by parts and by a Fourier expansion of the 
Kohn-Sham orbitals <^ a (x) m terms of plane waves u k {x). 
Thus, DFT practitioners have dealt with this form of 
nonlocality for a long time. This insight opens the av- 
enue for more complex (and hopefully more precise) den- 
sity functionals for nuclei. Recall that the very accurate 
mass formula by Duflo and Zuker |35| explicitly depends 
on the number of valence nucleons outside of closed shells. 
These are occupation numbers of shell- model orbitals, 
e.g., orbitals u k (x) from the harmonic oscillator or the 
Woods-Saxon Hamiltonian. Equation 14311 might permit 



us to explore and extend this very successful mass for- 
mula in the DFT setting. 

The main advantage of the functionals 1(12(1 and 1(39(1 
over the HFB approach consists of the conservation of 
particle number. This is very encouraging, as number 
projection techniques are numerically expensive and pose 
practical and theoretical difficulties [3(j. However, the 
Hamiltonian Q is not realistic, and a constant pairing 
interaction is a crude approximation for nuclei. Further- 
more, the involved sums over single-particle orbitals do 
not converge for infinitely large model spaces. To avoid 
divergences, one might introduce a cutoff (e.g., the Debye 
frequency for condensed matter problems), or utilize reg- 
ularization and renormalization procedures that relate 
the coupling g to physical observables |3H l38l l39l |40| . 
These methods not only cure divergences, they also sug- 
gest how to make the density functionals (|f 2(1 and 139(1 
more realistic. The density functionals of the pairing 
Hamiltonian are totally symmetric in the occupation 
numbers n k . This symmetry, of course, reflects the fact 
that the pairing interaction Q is a constant and equally 
scatters pairs between the single-particle orbitals. This 
suggests to modify the functional l|12(l through 

y/wfc(2 - n fc ) -> y^w fc \/n fc (2 - n k ) . (47) 

k k 

Here, uu k is an orbital-dependent weight or cutoff func- 
tion. The modified functional still exhibits the hallmarks 
of pairing, as they are described below Eq. (fT2"|) : appro- 
priately chosen weight functions w k also regularize the 
sums over single-particle orbitals, and should be capable 
of describing more realistic situations. 

VI. SUMMARY 

We constructed the density functional of the pairing 
Hamiltonian. This construction is possible in the strong 
coupling limit, and within the BCS approximation. In 
the strong coupling limit, the functional is a multiple 
product of local functionals that depend on the central 
moments of the occupation numbers. In the BCS approx- 
imation, the functional is arguably the simplest nonlocal 
and nonanalytical functional that causes an instability 
of the Fermi surface. Both functionals are based on oc- 
cupation numbers of single-particle orbitals that define 
the pairing Hamiltonian. These occupation numbers are 
themselves functionals of the Kohn-Sham orbitals that 
define the density. Our results suggest a way to include 
pairing with particle-number conservation into nuclear 
density-functional theory. 
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